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ABS TRACT: 


The response of a plane Poiseuille flow to disturbances of 
at «| various initial wavenumbers and amplitudes is investigated by 
. numerically integrating the equation of motion. It is shown that 
for very low amplitude disturbances the numerical integration scheme 
vi it yields results that are consistent with those predictable from 
a; linear theory. It is also shown that because of non-linear interactions 
a growing unstable disturbance excites higher wavenumber modes which 
have the same frequency, or phase velocity, as the primary mode. 
a For very low amplitude disturbances these spontaneously generated 
4 Ii higher wavenumber modes have a strong resemblance to certain modes 
rae computed from the linear Orr-Sommerfeld equation. 
7 


Brg ess § In general it is found that the disturbance is dominated for a 
| long time by the primary mode and that there is little alteration of 
the original parabolic mean velocity profile. There is evidence of 

’ the existence of an energy equilibrium state which is common to all 
Ly finite-amplitude disturbances despite their initial wavenumbers. 
This equilibrium energy level is roughly 3-5% of the energy in the 
et mean flow which is an order of magnitude higher than the equilibrium 
> value predicted by existing non-linear theories. 
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1. Introduction 


In this paper we present some results from a numerical investigation 
of the phenomena of stability, transition and turbulence in incompressible 
channel flows. Ideally, the goal of such an investigation would be to | 
numerically integrate the Navier-Stokes equations and thereby be able 
to follow the growth of some unstable disturbance which perturbs an 
initial laminar Poiseuille flow. It would then be possible to observe 
in detail how this growing disturbance eventually produces a fully 
developed, stationary turbulent flow. Presumably, if one were to repeat 
such a calculation for a large number of initial disturbances then the 
final turbulent flow could be described by an ensemble average over all 
the initial conditions. 

The results of such a computation would represent an exact (to 
within the limits of the numerical model) solution to the turbulence 
problem. Since turbulence is an inherently three-dimensional phenomenon, 
these computations would have to be carried out on a three-dimensional 
grid; unfortunately, the computational grid would have to be large enough 
to contain the large "eddies" representing the initial disturbance and 
still be able to resolve the smallest energy dissipation length scales 
of the final turbulence. Even at low Reynolds numbers the ratio of the 
larger length scale to the smallest scale is several orders of magnitude. 
The number of grid points is proportional to the cube of this ratio; 
hence, a complete three-dimensional solution is impractical, if not 
impossible, despite the large storage capacity of modern computers. 
(Emmons 1970) To overcome this impractically large storage requirement 


we have drastically simplified the problem by treating the flow on a 


two-dimensional basis. It is recognized that a two-dimensional treatment 
can not adequately represent true physical turbulence. Never-the-less 
the two-dimensional problem can be handled without further approximation 
and, as will be pointed out in the next section, two-dimensional 


solutions are of interest in their own right. 


2. Background 


During the past decade considerable effort was devoted to developing 
theories for the response of plane Poiseuille flows to finite-amplitude 
disturbances. The most significant contribution of these theories 
(in particular, those of Stuart and Watson (1960) and Reynolds and 
Potter (1967)) was in showing that the square of the amplitude ({a, |°) 
of an initially infinitesimal unstable disturbance is governed by an 
equation of the form 

d |A, \y 5 b 

ae Sy FG! he (1) 
Equation 1, which can be derived from the Navier-Stokes equations, is 
an approximation valid in a region in the (a-R, )-plane Which is close to 
the "neutral curve". Here R, denotes the Reynolds number and a the 
wave number of the disturbance. For disturbances with very small 
amplitudes the second term on the right of equation 1 must be neglible; 
hence, the square of the disturbance amplitude has an exponential 
growth and the constant, kK > is related to the exponential amplification 
factor of the linear theory (i.e., kK = -2 J B,,): Of particular interest 


is the case of k, > 0O for which the flow is unstable to small disturbances. 


if Ky, <0 then it is possible that the higher order terms will 


eventually balance the leading term and the amplitude growth is limited 


such that 
2 
d IA, | - 
at 
(1b) 
2 
var ait 
as 
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Under such conditions a "supercritical equilibrium state" is 
said to exist. On the other hand, if kK > O and k, > QO then no 
positive limiting value for la, IP is possible according to equation l. 
It is possible to evaluate Kk, based on Stuart's (1960) theory and 
the necessary calculations have been carried out by Pekeris and 
Shkoller (1967) and Reynolds and Potter (1967). These authors have 
found that under certain conditions k, <0O for some unstable disturbances. 
Both Pekeris and Shkoller and Reynolds and Potter have found that for 
regions near the lower branch of the "neutral curve" k, <0 while 
k, > O near the upper branch. Hence, one can conclude that disturbances 
whose wavelength corresponds to a point near the lower branch can grow 
according to the non-linear theory and reach a supercritical equilibrium 
state. On the other hand, for unstable disturbances whose wavelength 
corresponds to a point in the (a, R.) plane near the upper branch the 


non-linear theory does not seem to be applicable. At least the existing 


theories give us no information about the growth of such a disturbance 


beyond the linear range of amplitudes. 


Such "supercritical equilibrium states" are not observed experimentally 
in channel flows; instead the initial instabilities lead immediately 
to three dimensional turbulent flows. It seems, that no experimental 
verification of this aspect of finite-amplitude theory is possible. 

A computer simulation of a two dimensional channel flow is, perhaps, 
the best way of investigating the two dimensional response of a plane 
Poiseuille flow to an initial unstable disturbance. A computer 
simulation offers, in addition to its two dimensionality, the advantage 
of providing an "exact" solution to the equations of motion ("exact" 
to within the limits of a finite difference representation) without 
introducing any approximations concerning the relative magnitudes of 
the various modes of which the disturbance is composed. 

The present paper presents the results of a numerical investigation 
of the response of a plane Poiseuille flow to various initial disturbances. 
An attempt has been made to present the results in a manner that will 
facilitate comparison with the existing non-linear theory. The 
difficulty in making such comparisons is due to the fact that the 
theories mentioned concern the growth of an unstable disturbance with 
zero initial amplitude (at time t = -~) which grows at first according 
to linear theory and finally reaches an amplitude of such magnitude 
that the non-linear equations are appropriate. In a numerical simulation, 
with limited computation time available, the initial disturbance must 
have some finite amplitude. The growth rates predicted by the finite- 
amplitude theories are very small and the time required to follow 
the growth of a very small disturbance to the non-linear range consumes 


enormous amounts of computer time. For this reason the numerical 


investigator is tempted to start his calculations by using a disturbance 
whose amplitude is of such magnitude that the entire linear range can 
be bypassed. It will be pointed out in Section 6 that there seem to be 
Significant qualitative differences in the flow when large initial 
amplitudes are used as compared with those for which very small initial 
amplitudes are used in the calculations. 

Another aspect of the present investigation was to determine under 
what circumstances, if any, a two dimensional representation of the 
flow could be made to resemble in some respects three dimensional 
turbulent flows. 

It is clear that the "supercritical equilibrium states" in which 
the unstable disturbance mode remains dominant do not resemble 
turbulence in which the disturbance energy is distributed over a 
wide spectral range. Also in a "supercritical equilibrium state" 
the resultant mean velocity profile differs only slightly from the 
parabolic laminar profile while for turbulent flows the mean velocity 
profiles bear little resemblance to the original laminar one. It 
is recognized that there are significant differences between two 
dimensional and true three dimensional turbulence. This is particularly 
evident in the energy spectrum and mechanism for the transfer of 
energy and vorticity between the various spectral components (Lilly (1968) 
and Kraichnan (1967)). 

However, could it be possible that for some choice of initial 
disturbance mode shapes and amplitudes the structure of the resultant 


two dimensional flow could lie somewhere between the (relative) 


simplicity of a supercritical equilibrium state and the "complete 
chaos" characterizing a true turbulent flow? This aspect of the 


present paper will be discussed in Section 6. 


3. Basic Equations 


Consider the two dimensional flow of an incompressible fluid 


between two parallel planes. The governing equations of motion are 


2 


ae ay BC BY BCL 2 
at Oy Ox. Ox ay RY & (2a) 
and 
C= ¥ (2b) 


where € and Y¥ are the vorticity and stream function respectively with 


¢C and ¥ defined by 


a ee 7 = (3a) 
c= 2-2. (3b) 


The variables are assumed to be normalized on suitable characteristic 
lengths and velocities which are shown in Figure 1 along with other 


parameters of interest. The velocity components can be expressed as 


c 
Il 


Uy) #u 3; vyev 3s €2-=HrE 


(4) 


re 
[ 


= 3/2 (zy? -y) +¥' 


where the primed quantities which depend on x, y, t represent the 
instantaneous departure of the flow from the original laminar flow. 


Substituting equations 4 into equations 2 results in 
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ag" ag! ay' ay! 5 ' oy! 3 t 21 B) 
a a ets wR C (5a) 
c' = vy (5b) 


The conditions of zero slip at the wall and constant mass flow rate 


imply that the stream function, ¥', satisfy the boundary conditions 


y' = ayv'/’y =O at y=71. (6) 


By considering periodic solutions to equation 2 the boundary conditions 


in the streamwise direction are fixed by 


: i (x,y,t) y' (xf enL, Jy; t) (7a) 


4 


c’ (= Onn, ys 6) (Tb) 


UI 


. (yxsb) 


where the basic period is 2L. Equations 5a and 5b with the boundary 
conditions (6) and (7) are to be solved for ¥'(x,y,t) for a given 
initial value of ¥'(x,y,0). 

At this point it is convenient to define some parameters that are 
useful in describing the flow structure. The assumed periodicity 
implies homogeneity in the streamwise direction; hence, all averaging 
is done with respect to the streamwise coordinate, x. The average 
value of any quantity, q(x,y,t), is denoted by an overbar and obtained 


by integrating over one period of the basic wavelength. 


a 1: - 
a(yst)=s J a(x,y,t) ax (8) 
=f 


It should again be emphasized that in our representation the flow is 
arbitrarily divided into a laminar and a disturbance part. Hence, 


the total streamwise velocity component is given by 


u(e,y.t) - U(y) + a’ (xsy5t) (9) 
and the mean velocity is 
u(y,t) = U(y) + u” (y,t) (10) 


The difference between the mean and local values of any quantity is 
referred to as the turbulent part of that quantity and is denoted by 


a circumflex. Thus we have 


a Goyy t) = uilesy,t) - uly) =u" = 2 (11a) 
v (x,y,t) = v' (x,y,t) (11b) 
q (x75) = ¥' (xy¥,t) -¥ (y,t) (11c) 
a (arun) = ¢' ns t) = ee (yet) (11d) 
and 
A 
— = (1e) 
— _—ay' 
a a= 5 (1) 
Noy 
veo (11g) 
e.8. 
v= > ™ 0 (11h) 


The turbulent Kinetic energy is 


A re + 
B (yt) <3 [(- 2) + 7 (128) 


The mean value is 


a 1 Ly x 
f (y,t) = OL f E (x9 6) dx (12b) 
-L 


and the normalized sum of turbulent kinetic energy over a domain extending 
from the upper to the lower wall and for one period of the basic 


disturbance is 


A 1 1 en 
E(t) =3 f ay + 3 [ E (x,y,t) ax (12c) 
=i Ji 


Similarly, the kinetic energy in the mean flow is 


BE (y,t) =4 u(y) - aT (13a) 


and the normalized sum over the domain is 
1 
E(t)=3 | E (y,t) ay (13b ) 
-1 
The normalized energy in the initial laminar flow is 


B,=% f 30° ay = 600. (13c ) 


We can represent the periodic disturbance by means of a Fourier 


series as follows: 


Y Gets > £ot)e* (1a) 
n=-© 
¢' (x,y,t) = 5 g, (y,t) e nox (14 ) 
n==-o 
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with 


L 
al. inax 
rf (y,t) = OL, J * (x,y,t) dx (15a) 
= % 
fin = to 
1 - inax 
Bn (y,t) = oT i. C' (x,y,t)e dx (15b ) 
= % 
ea © a 


Where qa = n/L, and superscript * denotes a complex conjugate. From 5b 


2 
gs =f «(nay f. (15c) 


Notice that the mean value of the disturbance streamfunction and vorticity 


is simply 


¥' (yt) = £, (y,t) (16) 


Ul 


C' (y.t) = g (y,t) (17) 


The mean value of the turbulent kinetic energy is given by 


8 (y,t) = £ E (y,t) - E (y,t) (18) 


pe ©] 
where each of the components of the energy spectrum function, E(yst)s 


is given by 


2 ‘ 
E (yt) =3( |e) + @ : Ie, 1 ) (19) 


10 


The normalized total turbulent kinetic energy over the domain is 


R(t) = £ BL(t) - BA(t) (198) 


— Te) 


where 


b 
E(t) = 3 J E, (yst)ay (19b) 


4, The Numerical Model 


A rectangular grid of dimensions 2 x 2L is used with Mx N grid 
points located as shown in Figure 2. For the calculations presented 
herein 64 x 201 grid points were used. The indices k, j, and £ denote 
stations along x, y, and t respectively and 6x, Sy, 6&t are the 
corresponding intervals. The value of the disturbance streamfunction 
¥'(x,7,%) is then denoted by a - The use of the primes will be 
discarded from this point unless necessary for clarity. Equations 


5a and 5b are expressed in finite difference form as follows 


A411 p p p p 
So 7 “Kr, Ko “Kay Skeas ” K-27 
2 &t 26x J 2 6x 
(20a) 
L 1 2.4 
toes te * Ses 
4 2 6 2 ( p L p 
= =—— = + 
‘Ed “he oo VE “Gd “Reis 
(20b ) 
eo 4 4 4 ) 
oe he he” “teak 
by 


where A“ represents the modified form of the Laplacian differential 


operator of DuFort and Frankel (1958) defined as 


ae: 1_[/4 mf Atl 4 L 
/ en 
GK J WE Ska, 7 \Sxia * = or fa eae 
(21) 
of L ( Atl L- _ L ] 
a e 
aye | or a ek Sxya/ * Sk s-1 
and ae J is the total energy and mean square vorticity conservation 
form of the non-linear advection terms introduced by Arakawa (1966). 
Tis 
£1 {xy p - 
Ing * 5 Oe KT *y ge tee oe Se 
i: h h } 
§etcamb ster J 
Hs A Ce Akg AA AM g 


where 4. ant, 2, Gee iee ordinary first central differences. 


+ 
Equation 20a is solved explicity for os = at each time step from 
3 


£ h h-1 , 
the known values of "KJ ; GK J > and GK . The corresponding values 


of the Ye 3 are obtained from equation 20b. Equation 20b represents 
L 


a set of simultaneous algebraic equations to be solved for the by ay 
K=1,2,... M§ J=1,2,... N at each time step. Because of the large 
number of mesh points used in this investigation, conventional 
techniques are not suitable and an alternative method is used. We 


express the variables ne and Ye 5 in terms of their discrete Fourier 
3 b 


transforms as follows: 


L «i 
ws K 
ba? 2 tae (2ha) 
M 
h 4 = 8, K 
cKg7 % & 72 ? (2h ) 
With 
M 
<o,J0 = it 2 so : ° 
M 
A 1 £ ii 8n,x 
ei My OK & ? (24d) 
where 


Substituting equations 24 into equation 20b yields the following 


equation for the Fourier components of Yeo : 
bf 
% 2 
+ = 
foo” “yaa” “ayeet “ “y Ga,a (25) 
n = dae oM 
J me Doe eal 


where 


¥ 2 
a = 21 - (6y/&x)” (cos O42 - 1) | : 


The transforms and inverse transforms defined by equations 24 are 


computed by the fast Fourier algorithm of Cooley and Tukey (1965) while 


s 
gu 


equations 25 are tri-diagonal and are easily solved for the fs 


by Gauss elimination. 
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The assumed periodicity in x is automatically satisfied through 
equations 24, The boundary conditions at the walls (equation 6) are 


solved by setting ¥ - = ¥ 
K,1 


KN QO for all times and the no slip 
3 


condition is satisfied in the following way. We express Ys in 
3 


terms of ¥," through a Taylor series expansion about a point on the 


ea 
wall. 
¥ < Y ‘ + gy + A sy" +A. § 3 
ae et © hy Sone Ot ee 
£  .& e 3 
¥K,3 = Yt 2A, éy + 4A, éy + 8A. by~ + we. 
(26) 
4 ft 2 3 
Ku = Mea CA, 6y + oA,,, 6y + 27h. Cy oF ven 
oe 
The no slip condition (ay' /ay = 0) is satisfied if we require that 
A. = 0 which implies that (since ¥ - a 0) 
1 afal 
- 3» La 
2 7 2 %,3 7 9 Kh (272) 
and 
We me a ee (27) 
K,N-l * “K,N-2 9 K,N-3 


for all k and £. The Poisson equation (5b) is satisfied (to second 


order) at the walls by 


p _ i p p 

Kg 13/2 % 37 V9 YicsS (28a) 
£4 1 p p 

CK WN > a {3/2 YK Neo ~ 4/9 eat 285) 
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The length of the time increment, 6t, is fixed by the semi-empirical 


stability limit 


6t = f 6x evf\ (u + u). 6y + vo ox | (29) 


where the velocity components (U+u") | to . are the maximum absolute 
magnitudes of velocity selected from a large sample of grid points. 
Numerical experimentation has shown that for values of f < .6 consistent 
results are obtained. Values of f near unity may lead to catastrophic 
instability in this explicit technique. 

Since equation 20a is not “self starting" (it requires three time 


levels), the computation procedure starts with the forward difference 


representation 
4+1 L L+1 A+1 £+1 4+1 
“tg ~ ‘Kg . + {3 "Kt g 7 *K-1,3 aah ext g 7 &k-1,g 
6t ae 2 8x J 2 6x 


+1 7 a a 
Po + A GK 


K,J ” R, 
(30) 
p L p p 
1 "elgg ~ ey Ktl.J eK-1,3 
+2 2 6x we 2 8x 
4 1.2 a | 
oJ 
°K BO ies 


Equation 30 is also used periodically during the calculations in order 
to suppress any instabilities associated with the central time differencing 
used in equation 20a. 

The computation procedure can now be outlined step-by-step. 

1. An initial distribution of Ye 3 is chosen and the corresponding 


a are computed from equation 20b. 
3 


ID 


2. The time increment, &t, is computed from equation 29. 


3. The values of ie are computed for all the interior points 


3 
(J = 3,4,...N-2) by iteratively solving equation 30. 


4, The Fourier coefficients, ty z are computed from equation 
2hd, equation 25 ni solved for the f a and the corresponding 
values for the y z found from equation 2ha for the interior 
points (J = 3,4,...N-2). 


5. ‘The remaining grid points for Yo J ay pd = 1,2.N-1,N) 


are computed from equations ty Ob, 28a,b and 20b. 


6. For each time step the procedure in steps 3, 4, and 5 is 
repeated except the explicit central difference equation 20a 


is used in step 3 rather than equation 30. 


{. Periodically, say every 50 time steps, a new time increment, 
ét, is computed (step 2) and a single time step using the 
implicit forward difference equation 30 is used. 

For all of the computations presented in this paper a grid with 

dimensions M = 64, N = 201 was used which required a computation time 


of approximately 17 seconds per time step. 


5. Linear Calculations 


In an investigation such as this it is desirable to make some 
comparison between the results of our numerical computations and some 
known solution to the Navier Stokes equations. The imposed restraint 
of two dimensionality precludes any comparison with experimental data 
and, of course, there are no exact solutions to the non-linear equations 
of motion. The only alternative is to consider the limiting case of 


infinitesimal disturbances for which the results from linear theory 


16 


are available. For this reason, a rather thorough investigation of 
the linear solutions to equations 5a and 5b was undertaken. For small 
amplitude disturbances the non-linear terms in the vorticity transport 


equation are neglected and the equations of motion become 


oe + U “ = = = — . (31a) 
€ 
¢' = a i (31b) 


The general solution to equations 31 can be expressed in the form 


oo oo -i(nax - B am) 
Y'(x,y,t) = = y £49) e (32a) 
n=-© m=l 


co co -i(nax - Bm’) 
C'(x,y,t) = £ E Oy) ¢ (32b) 
=-o m=L 
with 
! 2 
.. i. (n a) Rs 
(32c ) 
a = 2n/r 


where } is the wavelength of the disturbance. 
Substituting equation 32 into 31 results in the well known Orr-_ 


Sommerfeld equation for $m) 1 Qs 


IV i 


B 
oa 2(na)° sf (ues g ting, ( a a2 \(5 


(33a) 
- (na)* 5) +U on | = 0 


ef 


with boundary conditions 
S = § = O at ye 2 Lx (33b ) 


solutions of 33a for one) for given a, n, and R. can be made to satisfy 
the homogeneous boundary conditions (33b) only for the discrete eigen- 
values, a 1smsoe, For the given a, n, and R. we are usually 
interested in determining if there exist corresponding Bam © which 
have negative imaginary parts. The corresponding disturbance, ome)? 
then grows exponentially with time and the flow is unstable to disturbances 
at the specified a, n, and Ro: 

For the present analysis we are interested in solutions to equation 
20a; hence, we must consider the finite difference representation of 


the linearized equation. For small disturbances equation 20a becomes 


41 bl p t t p 
7" SK _ , Ka “Ky Ses Kg 
26 : 2 x z 2 x (3h) 
2 >t 
+ yh 
/ e OK J 


Again we will consider solutions of the form (32); hence, we have 


~i(na(k-1) 6x - Bam & &t ) 


4 
ae = CS es - 
(2p) 
4 ~i(na(k-1) & - Cig * &t ) 
KI el nm © 


Substituting into 34 and 20b results in 


' . { iL 
Pram! °F) am = 196 [uj(@,) + (8) | ¥ R by° (85) 4 ae) oat ah _] 
e 
(362 ) 
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(8.) = — (a) 4 w ay (8,4) ] (36b ) 


nm Sy nm 

where 
w = 2((ay/ ex)” cos no ex) (37a) 

' 5 2 

wo = -2[1 + (=) (1 - cos na 6x) | (37b ) 
w' = sin na 6x/néx (37¢c) 

and 

; D sin B._ &t 
Bs 22 + 1/R, 6x | cos Bot +4 —= (374) 


i= /-1 


When expression 36b for (6,) is substituted into equation 36a there 
nm 


results a homogeneous set of linear algebraic equations for the (85) 
For given values of a, n, R.> 6x, Sy, and §t we wish to determine the 
eigenvalues and eigenvectors of the matrix of coefficients of the (85) 
In particular, we are interested in those eigenvalues Bs for which the 
imaginary part of een is negative and in the associated eigenvectors 

(35) The effects of 6x, Sy, and 6t on the linear solution are especially 
important for the present application. 

The eigenvalues and eigenvectors have been computed for a number of 
Reynolds numbers and for a wide range of a's using the QR algorithm 
(Wilkinson 1965 or Fair 1971). For more details on these calculations 
the reader is referred to O'Brien (1970) or to the report by Gawain and 


Clark (1971). In the latter report some of the eigenvalues Pa! and 
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eigenfunctions (ery) are tabulated for the case a = 1.0, R, = 6667 . 
In the limit as 6x, 6t ~ 0, the problem is the same as that solved by 
Thomas (1953) although the methods used are quite different from those 
employed by Thomas. For the case qa = 1.0, R= 6667, and éy = .01 (in 
Thomas' notation this corresponds to w= 1.0, R,= 10,000, éy = .01) we 
have solved equations 36 for n=1, 2, 3, and 4. Our solutions show 
that there is only one unstable eigenvalue and this unstable solution 
occurs for n=1. We shall order the eigenvalues, B am? so that for a 
given n, the m's are arranged according to the stability of the computed 
eigenvalues. Hence, for the case a = 1.0, R, = 6667, B,, is the unstable 
eigenvalue with the corresponding eigenfunction $1, (y). Thomas gives a 
value of B,, = .3653 - i .0055 and the corresponding mode shape, &4 Cy) 
is tabulated in his Table V. (Using our non-dimensionalization, B,, is 
exactly 1.5 times the value given by Thomas.) Our calculations yield 
the value Baa = .3593 - i .OO41 and the corresponding mode shape, %,5(y); 
is shown in Figure 3 for comparison with that of Thomas. 
Now consider the influence of the finite difference parameters 
&x, Sy, and §t on the solutions to 36. The effect of sy on the growth 
rate (ad B,4) is illustrated in Figure 4. As $y increases the flow 
becomes more and more stable and when 6y = .05, the phenomenon of linear 
instability has been completely obscured by the crude finite difference 
representation of the equations of motion. Indeed, calculations with 
$y = .05 for a wide range of Reynolds numbers and q's have shown no 
indication of linear instability. 
The effect of varying 6x was investigated at R, = 6667, a = 1.0, 


é6y = .O1 and &t < .02. Figure 5 shows the variation with 6x of the growth 
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rate (J B,,) and phase velocity (A, ,) for the unstable mode. For 
&x < .8 the changes in the disturbance profile (@,,(7)) where insignificant. 
Perhaps the most surprising discovery in this investigation of the 
linear behavior of the discrete representation of the equations of motion 
was the significance of the time interval, &t, on the disturbance growth 
rate (J B,,)- The effects of 6t on this parameter are shown in Figure 6. 
There are significant departures from the “exact” value of \) Boa for 
6t > .04 while the phase velocity and mode shape are not significantly 
changed. Hence, in addition to the usual stability constraint (29) on 
the time interval, &t, it is necessary that 6t < .O4 in the numerical 
integration scheme described in section 4. 
We can now summarize the results presented in the preceeding discussion. 
Increasing the y increment, Sy, tends to make the flow more stable and 
for coase grids (fy > .05) the phenomenon of linear instability is 
completely obscured. The linear solution is comparatively insensitive 
to changes in 6x, at least, until the x-increment is approximately the 
same length as the channel half-width. The growth rate of the disturbance 
is strongly influenced by the 6t increment (due to the DuFort-Frankel 
representation of the diffusion term) and for large 6t the growth rate 


is much higher than that predicted by the “exact” linear solution. 


6. Numerical Integration of the Vorticity Equation 


(a) Small Amplitude Disturbances 
We now return to the numerical integration of the vorticity 
equation using the techniques described in Section 4, As a check on 


the validity of the solution we first consider the limiting case of 


eal 


small amplitude disturbances and compare the results obtained with the 
known results from the linear solution as described in the preceeding 
section. To be more explicit, for a given Reynolds number and initial 
disturbance streamfunction, ¥' (x,y,0), are the results from the numerical 
integration consistent with the values of Band te) predicted by 
linear theory? 

A series of calculations were made at R, = 6667 and the initial 


disturbance given by 
' -i0x * +i ax 
¥' (x,y,0) = € {f,(y,0) e*™ + £,"(y,0) e** &*} (38) 


where ¢ is an amplitude factor. For e sufficiently small the non-linear 
terms in equation 5a should be neglible and the numerical integration 

should give results which are in agreement with the linear theory. As 

a first test we choose a= 1.0, R, = 6667, € = .005 and: the initial mode 
shape identical to the unstable mode computed from linear theory; i.e., 

f) (y,0) = $,,(y), where $,,(y) was shown in Figure 3. The growth 

rate of the disturbance kinetic energy, E (t), as defined by equation 12c 

is shown in Figure 7. The computed growth rate gives a value offs, = .0050 
which is in good agreement with the value given by Thomas. From equation 14 
we can determine the phase, @, (yt), of any of the spectral components 


' 
or ¥ Get) by 


@, (vst) = tan (J f (39) 


) 
n/pf 
The phase velocity, Cw or rate of propagation, of any of the spectral 


components is then 


o = + “h/ae (40) 
n nay 
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For these computations the grid length, 2L, was set equal to the wavelength, 
A, of the disturbance; so for the disturbance (primary) mode, n = 1, and 

in the linear range we can identify Cy with @ B,, from the linear theory. 
Figure 10 shows that Cy is, indeed, a constant and the slope of this 


phase versus time plot shows that c, = .363 which is in excellent 


li, 
agreement with Thomas' value of R B54 = .3653 and with our value of 
R B,, = 3593. 

Other short computations were carried out with, again, a = 1.0, 
a 6667, e = .005, but with the initial mode shape, f, (y,0), chosen to 
correspond to one of the stable modes computed from linear theory. These 
calculations were carried out by arbitrarily setting f, (y,0) equal to 

* 

the mode 7 $59: and $50 which are shown in Figure 9. The energy 
decay rates for these three cases are shown in Figure 10. From the energy 
decay rates and from the computed phase velocities (not. shown here) the 


numerical integration technique yields the following values for the 


appropriate eigenvalues 


Bae = 1.351 +i .138, Bog = 1.425 + i .207 


Boo * 1.456 + i .o76k. 


The linear matrix solution described in the preceeding section gave the 


values 


It 


By7 = 1.363 + i 136; Bog = 1.403 + i .1925 


P25 * 1.462 +4 .o7ue. 


* Since equation 33a is symmetrical in y it is possible to divide the 
solutions for @nm(y) into even and odd solutions. The boundary conditions 
are $ (+1) = $i n(tl) = 0, with S(O) = 6 (0) = 0 for even solutions and 
Bull = 4 (0) = « © for odd solutions. For our numerical model with N 
lateral Ba there are N-2 unknown (B35) v5 "s; hence, N-2 eigenvalues 
are computed for each n, a, R,. Of these N-2 eigenvalues, ay Coneeapoml 


to even solutions and N-3 correspond to odd solutions. 
2 
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A more critical test of the numerical integration is the ability to 
predict the correct mode shape, ODE for a disturbance of known wave— 
length and at a fixed R.° Suppose in equation 38 we choose a = 1.0 
but arbitrarily select the initial disturbance mode shape f, (y,0). The 
arbitrary function f, (y,0) can be expressed in terms of the eigenfunctions, 
ACRE of the Orr-Sommerfeld equation (33a) for the chosen a and R. 


(Schensted, 1961). i.e., 


f,(y,0) = 7" ey) (41a) 
with 
i 
A= J f(y,0) x,(y) dy (41 ) 
and 
y) = (25 - &) 3 
Xn cs 1m 


where 3, _(y) is the solution to the adjoint of equation 33a. However, 
the linear solution tells us that only one of the em in equation 41a is 
unstable; all other modes will decay with time. Hence, in the linear 
range only one mode, $4(y); grows and will eventually dominate all the 


other terms in the expansion (41a). Therefore, the mode shape 
(vst) +&,(y) as t >= 


where $4 (y) is as shown in Figure 3 for the case a= 1.0, R, = 6667. 
The initial shape chosen for f, (y,0) is shown in Figure 11. The 

Aw 
disturbance kinetic energy, E(t), is given in Figure 12 and the mode 


shapes, fi (yt), are shown for various times during the calculations. 


eu 


As expected, the disturbance energy at first rapidly decreases then 
grows according to linear theory while the computed mode shape is 
clearly tending towards the predicted shape of $4 y). 

From the results discussed in this section and from calculations 
at other wavenumbers, a, we have found that the numerical integration 
scheme consistently gives results that are in excellent agreement with 
the linear theory. The agreement between the results of the numerical 
model and the predicted results for small amplitude disturbances is 
gratifying and gives some confidence in the non-linear calculations 
now to be discussed. 

b) Finite Amplitude Disturbances 

Before discussing the results of the non-linear aspects of this 
investigation it is appropriate to give a very brief outline of the 
features of the existing non-linear theories which predict equations 
of the form of equation 1. The following discussion is from the 
paper by Pekeris and Shkoller (1967) which in turn aeiaine Stuart's 
(1960) theory for finite amplitude disturbances. If the Fourier 
representations for the disturbance streamfunction, ¥', and vorticity, 
C's; (equations 14 and 15) are substituted into the equation of motion 


(5a) the following equation for the modes results: 


og if t tf a | 
BR - (na)* g,) t+ina[(U-f, Je, - (U'-£,) fm] — (he) 


= ia om Oph pegs ac 


where 


i (na) fm (15c ) 


oy) 


and 


Ay = {(n 8) (f, En-s ~ &s f-s) 


(42a ) 


. (nts)(£," Ents 7 g. a) 


if the parameter, & = R| B,, |, is small (implying proximity to the 
neutral curve) then Stuart showed that in an approximation in which 
in (14a) 


terms of order - 3/2 are retained, only the f tf)» and f 


O* 2 
are significant. Now expand f, (yt) in terms of the eigenfunctions, 
$ of the Orr-Sommerfeld equation as was done in equation (41a), 


so that 


fi(y,t)=¢ EF A(t) 5 (43) 
m=1 


Stuart's asymptotic analysis for small e. shows that only the first 


term in (43) is important, i.e., 
fi(y,t) # e, A(t) ,(y) +0 (e°) (Iu) 


and that the terms f and f, are of smaller order so that 


, 2 
e(y,t) = ela (t)| Gly) + 0(e,*) (45) 


e.(yst) = ¢ aM(t) Gy) Of) (46) 


The functions G. and G, are determined from Pekeris and Shkoller's 
equations (22), (23), and (24). By substituting equations (44), (45), 
and (46) into (42) one can obtain an expression of the form (1) for 


ed 
the amplitude parameter, JA, (t) | 5 Sins Bing 
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aay 
= ra +k, ran (1) 
The important point from our viewpoint is that according to the above 
analysis, the disturbance energy should be contained primarily in the 
f(y,t) mode and that the growth of this energy should behave according 
to equation (1). Assuming that the disturbance streamfunction can be 


adequately represented by f, then we have from (44) 


1 


-i ax 


¥'(x,y,t) = ¢ (A,(t) &,(y)e™ %* + ap (t) 45 (ye °*) (47) 


The disturbance kinetic energy as defined by equation (12c) (since we 


are neglecting the contribution from f, the "disturbance" and "turbulent" 


6) 


kinetic energies are the same) becomes 


E(t) = 6° T a(t) | (48) 
where 


I = f fies t 7's is,.1 } dy (49) 


It was pointed out in part 2 of this paper that the constant, > ax 
(1) is related to the amplification factor of linear theory (k,= -2\ Pag) 
and if ko <0O then the growth of an unstable disturbance is limited 


so that 
als 1/k, (50) 


Pekeris and Shkoller found that these conditions are satisfied for a 
and R. inside the lower portion of the neutral curve. For example, 


at R, = 6667 and a = .90 their calculations give 
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2 
4,1 71/9 


Our solutions to the linear eigenvalue problem at ‘ie 6667 and 


a = .875 gives BL = -.00394 and I = 1.848. Hence, 
a lb 
E(t) > (1,848)(.00394)/9.0 = 8.08 x 10° 
Figure 13 shows the energy growth for modes Eo» Ea: E,> and B. for 


a run with the initial disturbance given by (38) with a = .875, 
R, = 6667, € = .005, and f, (y,0) = $,,(y) where $,, in the unstable 
eigenfunction for a = .875. The disturbance is, indeed dominated by 
B, (t) and the mode shape, f, (y,%), normalized so that f, (0,t) = 1.0, 
never departed from the initial shape given by $,,(y)- A longer 
calculation was made using the same initial conditions except in this 
case w=1.0. Again the energy remained predominately in the n =1 
mode and again the mode shape for f, (y,t) never departed from the initial 
$4 (y) shape indicating that (44) and (47) do indeed adequately represent 
the flow as far as the overall disturbance energy is concerned. 

Although this latter calculation required 34.5 hours of computer 
time, we can see that the energy growth is still within the linear 
range. This demonstrates that, as was pointed out in part 2, it is 
impractical using the present techniques to follow the disturbance 
through the linear and into the non-linear range of amplitudes. For 
this reason we have been unable to establish whether or not the 
"supercritical" equilibrium states as predicted by (1) and (50) do, 
in fact, exist for our exact treatment. Our calculations, as shown 


in Figures 13 and 14, do indicate, however, that there is no significant 
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difference between cases for which a = .875 or a = 1.0 in equation (38). 
The calculation of Pekeris and Shkoller, however, show that “super- 
critical" equilibrium states are possible for a = .875 but not for 
a= 1.0. All our calculations for both low and high amplitude disturbances 
have shown no strong dependence on a as far as the final result is 
concerned. 

Although this last run was entirely within the linear energy 
growth range we will see that there were very interesting non-linear 
features to the calculations. The appearance of energy in modes 
9? Ep» and E, as shown in Figure 14 is, of course, due to the 


presence of the non-linear terms in equation (5a). In identical 


E., EB 
calculations for which the non-linear terms in (5a) were deliberately 
omitted the primary (n=l) mode grew according to linear theory while 
the energy in the other modes remained constant or decreased and the 
energy levels (due to numerical “noise") were twelve to fifteen orders 
of magnitude lower than that of the primary mode. Figure 14 shows 
that the energy in the other modes (n = 0,2,3) grows very rapidly at 


first and then seems to follow an exponential growth rate, i.e., 


eB 


E(t) ~ 2 n = 0,1,2,3 


for large t. From the data shown in Figure 14 the values for B 


are By = Olt, B, = .0050, B, = .0080, and B. = .0117. 
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The most interesting feature is the shape of the modes, f,» f3. 
with respect to y. These shapes are shown in Figure 15 where the 
modes are normalized so that f (0,t) = 1.0 for even modes and for 
odd modes the maximum real and imaginary parts of f(y) are unity. 
The phase velocities for these modes averaged between t = 65 and 


t = 71 were c, = .362 and c.. = .363. These phase velocities are the 


2 
Same as the phase velocity of the primary mode predictable from linear 
theory (i.e. 6,, = .363). 

The striking feature is the strong resemblence between these modes 
and the modes calculated from the Orr-Sommerfeld equation. Indeed, it 
appears that the non-linear solution has spontaneously generated modes 
that can be predicted from the linear theory even though the linear 
solution indicates that these are stable modes and would decay with 
time. This means that if we represent each mode by an expression of 


the form of (43) i.e., 


f (y,t) = 3 A(t) SW) (51) 


*From equation (15c) we see that g, has the same symmetry or anti- 
symmetry as does f,3 i.e. if f,(y) is an even function of y then gy 
is also even. Then from equation (42) we can establish that if 
¥'(x,y,0) is given by (38) and f,(y,0) is even in y then fo(y,t) is 
odd, fo(y,t) is odd, f3(y,t) is even, fy(y,t) is odd, etc. Likewise 
if f1(y,0) is odd then fo(y,t) is odd, fo(y,t) even, f3(y,t) odd, 
fi(y,t) is even, etc. For example, let fi(y,0) be even. Then since- 
f_(y,0) #0 ‘fp * 1 


= odd function of y. 
At the next time step, t = &t 


g(y, ét ) - P&2/at) st = odd function of y 
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then one mode seems to dominate the expansion. It is interesting to 
notice that in each case the dominant mode is one for which there exist 
a step gradient near the wall. Most of the stable modes that we have 
computed from linear theory are similar to the ones shown in Figure 9 
which have a detailed structure near the center of the channel rather 
than near the wall. 

A series of higher energy runs with various initial disturbance 
amplitudes was made at R, = 6667 using as initial conditions a 
disturbance with wave number qa = 1.0 and that mode shape given by Thomas 
(see Figure 3). The total energy growth (or decay) is shown in Figure 16 
for several different initial energy levels. A further energy breakdown 
is given in Figure 17 for the cases e¢ = (.05)/ 5. For this case the 
total energy remains almost constant while there is a continual 
exchange of energy between the mean flow and the turbulence. The 
turbulent kinetic energy oscillates about a value of E(t) = .03 which 
is roughly ten times the equilibrium value predicted from equation 50. 
From Figure l4a it seems that there is a tendency towards some common 
energy equilibrium state for the different initial conditions. A more 
detailed inspection of the flow structure shows that there are significant 
differences between the three flows. Starting with the high energy run 
we summarize the structure of the flow by presenting plots of the mean 
velocity profile, Reynolds stresses, and energy spectrum at various 
times throughout the run (Figures 18, 19, and 20). Also, in Figure 21 
the shape of the primary (n = 1) mode at the end of the calculations 
is shown for each of the three energy levels. For the high energy run 


(although the time parameter is comparatively small) there are drastic 
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alterations of the mean velocity profile and primary mode shape and 
considerable energy transfer from the primary mode to the higher 
harmonics. For the two lower energy runs there is little alteration 

of the parabolic velocity profile and only small amounts of energy 
transfer to higher wave numbers. For the lowest initial energy 

level (e = .05 / 2) the primary mode shape remains essentially unchanged 
throughout the entire calculation. More details of this low energy 

run are shown in Figure 22 where the growth of the total turbulent 

energy in the first three spectral components (equation 19b) is 

shown. The phase velocities (equation 40) were computed near the 


end of the run and are summarized below: 


c, = 492, e4) = MOH, CG. = 473 
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The measured C_ were independent of y and all are roughly 4O% higher 
than the wave velocity (Ap,,) of the primary unstable mode predicted 
by linear theory. This low energy run has some qualitative similarities 
to the predictions of the non-linear theory which was outlined at the 
beginning of this section. That is, there appears to be an asymptotic 
energy level, the primary mode shape is only slightly altered and the 
energy in the primary mode remains dominant. 

In order to investigate the influence of initial wavenumber on 
the flow two runs were made at R, = 6667 with a = .875 (near the 
lower branch of the neutral curve) and with a = 1.05 (near the upper 
branch). For both cases the initial mode shape was determined from 
the eigenvalue solution to the linear equations and the initial 


energy was chosen to be slightly lower than that of the run (e = .05/) 5; 


Figure 16) which seemed to remain in an energy equilibrium state. The 
energy variation for these cases is shown in Figure 23. There is no 
obvious qualitative difference between the three wavenumbers; the initial 
high energy level, however, apparently causes oscillations of the 
disturbance kinetic energy and nakes any direct comparison with non- 
linear theory rather difficult. 

Another important question which arises in the consideration of 
non-linear aspects in the stability of plane Poiseuille flows is whether 
or not a flow which is stable to infinitesimal disturbances might be 
unstable to disturbances of some finite amplitude. Such instabilities 
are generally referred to as “subcritical” instabilities. To investigate 
this aspect of the problem several runs were made using as initial 
conditions a disturbance that is known to be stable in he linear 
range. The case chosen was for a = .78, i 6667 for which B,, = .2hkg 
+ 4 .00124 for the least stable eigenfunction determined from the linear 
solution described in Section 5. Again the initial disturbance was 
given by equation 38 with the function f, (y,0) = $4 (y) determined 
from the linear solution. For a small amplitude disturbance (e = .005) 
the disturbance energy decreased with time according to linear theory 
while for a larger amplitude (e = .0635). The energy variation was 
as shown in Figure 24. The familiar oscillation of the turbulent 
kinetic energy which is associated with an initially high disturbance 
energy is present in this case; while the total kinetic energy appears 
to remain approximately constant after an initial decrease. 

We now address ourselves to the question posed at the end of 


Section 2 concerning the character of a two dimensional "turbulence" 
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in channel flows. For all the runs discussed above the computational 
grid length (2L) was equal to the wavelength of the initial unstable 
disturbance (i.e. 2L = 2n/a). In any finite difference representation 
of the equations of motion there are only discrete wavelengths 
available to be excited by the non-linear interactions. Recall from 


equation 24a that the streamfunction is expressed as 


M -1¢, k 
¥ = > f e : ela 
KJ eT nod (2ha ) 
and since "x is real 
R Fad = ty-nt2 J 


J Aad . “J Fnen+2,J 


Lesh r is even and = t is odd over the interval n = 1,2,...M/2, 
n n 


J J 
3 2 
...M. Therefore, there are M/2 distinct Fourier modes "available" to 
the solution of the equations of motion. When the grid length equals 
the wavelength of the initial disturbance there is only one lower mode 


available to be excited; namely, the zero mode, f It is known 


0,J 
that two-dimensionality constrains the transfer of energy between 
the modes of which the disturbance is composed. For an isotropic 
two-dimensional turbulence there can be no systematic transfer of 
energy from lower to higher wave numbers without a corresponding 
transfer to still lower wavenumbers (Lilly (1968), Kraichnan (1967)). 
Hence, one may suspect that our discrete representation in which 


there are no available modes at lower wave numbers could force the 


disturbance energy to remain in the primary initial mode. The first 
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alternative to this dilemma would seem to be to increase the grid 
length, 2eL, thereby allowing modes with both higher and lower wave 
numbers to fit into the available grid length. However, if the 
initial disturbance energy is contained in, say, wave number n then 
the non-linear interaction of this mode with itself will initially 
excite wave numbers O and 2n; later modes 0, 2n, 3n, Un, etc. ‘The 
modes that can be excited will be separated by a spacing of n in 
wave number space, resulting in a spectrum with "holes" or regions 
in which no energy can appear. Clearly, this is not a realistic 
case. The alternative to this situation is to initially perturb the 
flow with a disturbance in which at least two contiguous (in the 
discrete wave number space) modes are present. For example, if the 
initial disturbance contains modes n and n+l then the non-linear 
interaction will, at first, excite modes 0, l, en, entl, ent2; these 
in turn will excite modes 0, 1, n-l, 2n-1, 4n, etc. until all possible 
wave numbers have been excited. 

It seems reasonable to require that the wavelengths of these two 
initial modes should both fall within the unstable region in the 
(a-R,) plane. If we set the grid length so that 2L=16n (or a = 1/8) 
then modes with wavenumbers 7a and 8q@ will both fall within the 


unstable region. Hence 


¥'(x,y,0) = 6 {tne 870% 2 — m oe 4 ian 


be 


A lengthy numerical integration of the vorticity equations using 
the initial conditions just described is summarized in Figures 25 
and 26 which show the turbulent kinetic energy time variation and 
the energy spectrum at various times during the run. For this run 
the turbulent kinetic energy initially grows and then appears to 
vary randomly about a value of E(t) = ,019. The energy spectrum is 
quite different from the other cases presented in that there is a 
gradual transfer of significant amounts of energy to other wavenumbers. 
Although it is not shown here, the mean velocity profile again differed 
only slightly from the parabolic distribution. Although, there 
are significant differences between this run and the others presented, 
all the examples that have been discussed are qualitatively similar 
in that most of the energy remains in modes whose wavelength 


approximately equals that of the primary (unstable) mode. 


7. Conclusion 


Our numerical technique for integrating the vorticity equation 
has been shown to be stable and predicts results that are consistent 
with linear solutions to the vorticity equation. 

The influence of the finite difference representation on the 
linear eigenvalue solution to the equations of motion was investigated 
and some interesting conclusions resulted. The size of the 6y 
increment significantly influences the stability of the equations 
and for grids with Sy > .05 (approximately) the equations are stable 
to all disturbance. Because of the DuFort-Frankel representation 
of the diffusion term, the growth rate of an unstable disturbance is 


strongly influenced by the 6&t increment. The growth rate increases 
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rapidly with 6t even though the time step increment may still be 
small enough to satisfy ordinary numerical stability requirements. 
The 6x increment has the smallest effect on the linear solution 
and seems to be significant only when 6x is approximately equal to 
the channel half-width. 

Because of the long computing time required we were unable to 
follow the growth of an unstable disturbance through the linear 
and into the non-linear range of disturbance amplitudes. This 
precluded any direct comparison of our results with those predicted 
by the non-linear theories of Stuart and Watson and Reynolds and Potter. 
We can, however, make some general comparisons. We find qualitative 
agreement in the sense that the unstable primary mode dominates the 
disturbance and there is little distortion of the mean flow. The 
calculations presented indicate the existence of a disturbance energy 
equilibrium state which seems to be common (for a given Reynolds 
number) to all unstable disturbances despite their initial energy 
or wavenumber. Indeed, we have found no outstanding differences 
in the behaviour of disturbances with varying initial wave numbers 
(at a constant R,) which is not in agreement with the non-linear 
theories. The distortion of the primary mode shape is significant 
when the initial energy is high enough to be well beyond the linear 
range but for lower initial energies the distortion of the mode 
shape is less important. 

The computed flows which have been discussed in this paper 
show little resemblance to actual turbulent shear flows. In a true 


three dimensional flow one would expect an initial unstable disturbance 
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with a length scale approximately 2n times the channel half width to 
grow in amplitude and through non-linear interactions the energy of 
the disturbance would be transferrred to higher wavenumbers with length 
scales comparable to the channel half width. The corresponding 
situation for a two dimensional turbulent channel flow is not clear. 
However, for two dimensional homogeneous turbulence the energy 
would cascade to the lower wavenumbers while mean square vorticity 
is transferred to higher wavenumbers when the turbulent energy is 
continually fed into the flow in a narrow band of spectral components. 
Neither of these processes was observed for the present calculatims. 
One can speculate that our calculated "quasi-equilibrium" states 
in which most of the energy remains in one or two spectral components 
with little alteration of the parabolic velocity profile represent 
two dimensional solutions to the equations of motion which are, in 
some sense, highly unstable flows and, for this reason, are not 


observed experimentally. 
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Figure 7. Energy Growth for a Low Amplitude Run. 
R, = 6667, a = 1.0 
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Figure 11. Primary Mode Shape, f,, for a Low Energy 
(Linear) Run at Times t=0 and t268, 
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Figure 19b. Reynold's Stresses for the Run 
e= .05, a = 1.0, R, = 6667 
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Figure 20a. Mean Velocity Profiles for the Run 
t= .05/,/2, aw = 1.6, R, = 6667 


67 


| 
| 
| 
| 
| 
| 
y | 
1 
0 | 
| 
| 
| 
| 
| 
—1.0 
-.0003 O AA .0003 
UV 
T=111.8 
1.0 
| 
| 
| 
| 
| 
| 
y | 
| 
0 
| 
| 
| 
| 
| 
| 
| 
| 
950003 O .0003 


Figure 20b. Reynold's Stresses for the Run e = .05/,/2 
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